/*!
 *  \file       PoissonRandomGenerator.cpp
 *  \brief      Implementazione della classe \c PoissonRandomGenerator
 *  \author     Davide Scola <S2299453@studenti.unige.it>
 *  \author	Minetti Alberto <S3064777@studenti.unige.it>
 *  \version    1.0
 */
#include <cmath>
#include "IEntropyGenerator.hpp"
#include "IRandomGenerator.hpp"
#include "SimpleRandomGenerator.hpp"
#include "NormalizedRandomGenerator.hpp"
#include "PoissonRandomGenerator.hpp"

using namespace std;


/*!
 *  \brief      Costruttore di copia di \c PoissonRandomGenerator
 *  \param[in]	other	generatore da copiare
 */
template < typename tType >
PoissonRandomGenerator< tType >::PoissonRandomGenerator( PoissonRandomGenerator< tType > const& other )
  : NormalizedRandomGenerator< tType >( static_cast< NormalizedRandomGenerator< tType > >( other ) ),
    mean_( (&other)->mean_ )
{
}

/*!
 *  \brief      Costruttore di \c PoissonRandomGenerator
 *  \param[in]	mean		media
 *  \param[in]	entropy		generatore di entropia da utilizzare
 */
template < typename tType >
PoissonRandomGenerator< tType >::PoissonRandomGenerator( tType mean, IEntropyGenerator & entropy )
  : NormalizedRandomGenerator< tType >( entropy ),
    mean_( mean )
{
}


/*!
 *  \brief      Operatore di assegnazione della classe \c PoissonRandomGenerator
 */
template < typename tType > PoissonRandomGenerator< tType >&
PoissonRandomGenerator< tType >::operator=( const PoissonRandomGenerator& other )
{
   if (this == &other) {
       return *this;
   }
   
   this->entropy_ = (&other)->entropy_;
   this->mean_ = (&other)->mean_;
   
   return *this;
}


/*!
 *  \brief      Funzione ausiliaria per la distribuzione poissoniana
 *  \warning	Funzione solo adattata, presa da Aulaweb
 */
template < typename tType > tType
PoissonRandomGenerator< tType >::gammln( tType xx )
{
  tType x, tmp, ser;
  static double cof[6] = {
    76.18009173, -86.50532033, 24.01409822, -1.231739516, 0.120858003e-2, -0.536382e-5
  };
  
  ser = tType( 1 );
  x = xx - tType( 1 );
  tmp = x + tType( 5.5 );
  tmp -= ( x + tType( 0.5 ) ) * tType( log( tmp ) );
  
  for( unsigned int j = 0; j <= 5; ++j ) {
    x += tType( 1 );
    ser += ( (tType) cof[j] ) / x;
  }
  
  return -tmp + log( tType( 2.50662827465 ) * ser );
}

/*!
 *  \brief      Ritorna un nuovo numero pseudo-casuale con distribuzione poissoniana
 *  \warning	Funzione solo adattata, presa da Aulaweb
 */
template < typename tType > tType
PoissonRandomGenerator< tType >::poidev( tType xm )
{
  tType em, t, y;
  static tType sq, alxm, g, oldm = tType( -1 );
  
  if( xm < tType( 12 ) ) {
    if (xm != oldm) {
      oldm = xm;
      g = tType( exp( -xm ) );
    }
    
    em = tType( -1 );
    t = tType( 1 );
    
    do {
      em += tType( 1 );
      t *= NormalizedRandomGenerator< tType >::next( );
    } while( t > g );
  }
  
  else {
    if( xm != oldm ) {
      oldm = xm;
      sq = tType( sqrt( tType( 2 ) * xm ) );
      alxm = tType( log( xm ) );
      g = xm * alxm - this->gammln( xm + tType( 1 ) );
    }
    
    do {
      do {
	y = tType( tan( M_PI * NormalizedRandomGenerator< tType >::next( ) ) );
	em = sq * y + xm;
      } while( em < tType( 0 ) );
      
      em = tType( floor( em ) );
      t = tType( 0.9 ) * ( tType( 1 ) + y * y ) * tType( exp( em * alxm - this->gammln( em + tType( 1 ) ) - g ) );
    } while( NormalizedRandomGenerator< tType >::next( ) > t );
  }
  
  return em;
}


/*!
 *  \brief      Ritorna un nuovo numero pseudo-casuale con distribuzione poissoniana
 */
template < typename tType > tType
PoissonRandomGenerator< tType >::next( )
{
  return this->poidev( this->mean_ );
}
